10 Area Data I

Area Data I

NOTE: You can download the source files for this book from here. The source files are in the format of R Notebooks. Notebooks are pretty neat, because the allow you execute code within the notebook, so that you can work interactively with the notes.

In last few practices/sessions, you learned about spatial point patterns. The next few sessions will concentrate on area data.

If you wish to work interactively with this chapter you will need the following:

• An R markdown notebook version of this document (the source file).

• A package called geog4ga3.

Learning objectives

In this practice, you will learn:

1. A formal definition of area data.
2. Processes and area data.
3. Visualizing area data: Choropleth maps.
4. Visualizing area data: Cartograms.

Suggested reading

O’Sullivan D and Unwin D (2010) Geographic Information Analysis, 2nd Edition, Chapter 7. John Wiley & Sons: New Jersey.

Preliminaries

As usual, it is good practice to clear the working space to make sure that you do not have extraneous items there when you begin your work. The command in R to clear the workspace is rm (for “remove”), followed by a list of items to be removed. To clear the workspace from all objects, do the following:

rm(list = ls())

Note that ls() lists all objects currently on the worspace.

Load the libraries you will use in this activity:

library(tidyverse)
#library(rgdal)
#library(broom)
library(plotly)
library(cartogram)
library(gridExtra)
library(geog4ga3)

Read the data used in this chapter.

data("Hamilton_CT")

The dataframe includes the spatial information for the census tracts in the Hamilton Census Metropolitan Area in Canada (as polygons), and a series of demographic variables from the 2011 Census of Canada.

You can quickly verify the contents of the dataframe by means of summary:

summary(Hamilton_CT)

       ID               AREA             TRACT             POPULATION     POP_DENSITY       
 Min.   : 919807   Min.   :  0.3154   Length:188         Min.   :    5   Min.   :    2.591  
 1st Qu.: 927964   1st Qu.:  0.8552   Class :character   1st Qu.: 2639   1st Qu.: 1438.007  
 Median : 948130   Median :  1.4157   Mode  :character   Median : 3595   Median : 2689.737  
 Mean   : 948710   Mean   :  7.4578                      Mean   : 3835   Mean   : 2853.078  
 3rd Qu.: 959722   3rd Qu.:  2.7775                      3rd Qu.: 4692   3rd Qu.: 3783.889  
 Max.   :1115750   Max.   :138.4466                      Max.   :11675   Max.   :14234.286  
  AGE_LESS_20      AGE_20_TO_24    AGE_25_TO_29    AGE_30_TO_34     AGE_35_TO_39   
 Min.   :   0.0   Min.   :  0.0   Min.   :  0.0   Min.   :   0.0   Min.   :   0.0  
 1st Qu.: 528.8   1st Qu.:168.8   1st Qu.:135.0   1st Qu.: 135.0   1st Qu.: 145.0  
 Median : 750.0   Median :225.0   Median :215.0   Median : 195.0   Median : 200.0  
 Mean   : 899.3   Mean   :253.9   Mean   :232.8   Mean   : 228.2   Mean   : 239.6  
 3rd Qu.:1110.0   3rd Qu.:311.2   3rd Qu.:296.2   3rd Qu.: 281.2   3rd Qu.: 280.0  
 Max.   :3285.0   Max.   :835.0   Max.   :915.0   Max.   :1320.0   Max.   :1200.0  
  AGE_40_TO_44     AGE_45_TO_49    AGE_50_TO_54    AGE_55_TO_59    AGE_60_TO_64  AGE_65_TO_69  
 Min.   :   0.0   Min.   :  0.0   Min.   :  0.0   Min.   :  0.0   Min.   :  0   Min.   :  0.0  
 1st Qu.: 170.0   1st Qu.:203.8   1st Qu.:203.8   1st Qu.:175.0   1st Qu.:140   1st Qu.:115.0  
 Median : 230.0   Median :282.5   Median :280.0   Median :240.0   Median :220   Median :157.5  
 Mean   : 268.7   Mean   :310.6   Mean   :300.3   Mean   :257.7   Mean   :229   Mean   :174.2  
 3rd Qu.: 325.0   3rd Qu.:385.0   3rd Qu.:375.0   3rd Qu.:325.0   3rd Qu.:295   3rd Qu.:221.2  
 Max.   :1105.0   Max.   :880.0   Max.   :740.0   Max.   :625.0   Max.   :540   Max.   :625.0  
  AGE_70_TO_74    AGE_75_TO_79     AGE_80_TO_84     AGE_MORE_85              geometry  
 Min.   :  0.0   Min.   :  0.00   Min.   :  0.00   Min.   :  0.00   POLYGON      :188  
 1st Qu.: 90.0   1st Qu.: 68.75   1st Qu.: 50.00   1st Qu.: 35.00   epsg:NA      :  0  
 Median :130.0   Median :100.00   Median : 77.50   Median : 70.00   +proj=long...:  0  
 Mean   :139.7   Mean   :118.32   Mean   : 95.05   Mean   : 87.71                      
 3rd Qu.:180.0   3rd Qu.:160.00   3rd Qu.:120.00   3rd Qu.:105.00                      
 Max.   :540.0   Max.   :575.00   Max.   :420.00   Max.   :400.00                      

Area data

Every phenomena can be measured at a location (ask yourself, what exists outside of space?).

In point pattern analysis, the unit of support is the point, and the source of randomness is the location itself. Many other forms of data are also collected at points. For instance, when the census collects information on population, at its most basic, the information can be georeferenced to an address, that is, a point.

In numerous applications, however, data are not reported at their fundamental unit of support, but rather are aggregated to some other geometry, for instance an area. This is done for several reasons, including the privacy and confidentiality of the data. Instead of reporting individual-level information, the information is reported for zoning systems that often are devised without consideration to any underlying social, natural, or economic processes.

Census data, for instance, is reported at different levels of geography. In Canada, the smallest publicly available geography is called a Dissemination Area or DA. A DA in Canada contains a population between 400 and 700 persons. Thus, instead of reporting that one person (or more) are located at a point (i.e., an address), the census reports the population for the DA. Other data are aggregated in similar ways (income, residential status, etc.)

At the highest level of aggregation, national level statistics are reported, for instance Gross Domestic Product, or GDP. Economic production is not evenly distributed across space; however, the national GDP does not distinguish regional variations in this process.

Ideally, a data analyst would work with data in its most fundamental support. This is not alway possible, and therefore many techniques have been developed to work with data that have been agregated to zones.

When working with areas, it is less practical to identify the area with the coordinates (as we did with points). After all, areas will be composed of lines and reporting all the relevant coordinates is impractical. Sometimes the geometric centroids of the areas are used instead.

More commonly, areas are assigned an index or unique identifier, so that a region will typically consist of a set of \(n\) areas as follows: \[ R = A_1 \cup A_2 \cup A_3 \cup ...\cup A_n. \]

The above is read as “the Region R is the union of Areas 1 to n”.

Regions can have a set of \(k\) attributes or variables associated with them, for instance: \[ \textbf{X}_i=[x_{i1}, x_{i2}, x_{i3},...,x_{ik}] \]

These attributes will typically be counts (e.g., number of people in a DA), or some summary measure of the underlying data (e.g., mean commute time).

Processes and area data

Imagine that data on income by household were collected as follows:

df <- data.frame(x = c(0.3, 0.4, 0.5, 0.6, 0.7), y = c(0.1, 0.4, 0.2, 0.5, 0.3), Income = c(30000, 30000, 100000, 100000, 100000))

Households are geocoded as points with coordinates x and y, whereas income is in dollars.

Plot the income as points (hover over the points to see the attributes):

p <- ggplot(data = df, aes(x = x, y = y, color = Income)) + 
  geom_point(shape = 17, size = 5) +
  coord_fixed()
ggplotly(p)

The underlying process is one of income sorting, with lower incomes to the west, and higher incomes to the east. This could be due to a geographical feature of the landscape (for instance, an escarpment), or the distribution of the housing stock (with a neighborhood that has more expensive houses). These are examples of a variable that responds to a common environmental factor. As an alternative, people may display a preference towards being near others that are similar to them (this is called homophily). When this happens, the variable responds to itself in space.

The quality of similarity or disimilarity between neighboring observations of the same variable in space is called spatial autocorrelation. You will learn more about this later on.

Another reason why variables reported for areas could display similarities in space is as an consequence of the zoning system.

Suppose for a moment that the data above can only be reported at the zonal level, perhaps because of privacy and confidentiality concerns. Thanks to the great talent of the designers of the zoning system (or a felicitous coincidence!), the zoning system is such that it is consistent with the underlying process of sorting. The zones, therefore, are as follows:

zones1 <- data.frame(x1=c(0.2, 0.45), x2=c(0.45, 0.80), y1=c(0.0, 0.0), y2=c(0.6, 0.6), Zone_ID = c('1','2'))

If you add these zones to the plot:

p <- ggplot() + 
  geom_rect(data = zones1, mapping = aes(xmin = x1, xmax = x2, ymin = y1, ymax = y2, fill = Zone_ID), alpha = 0.3) + 
  geom_point(data = df, aes(x = x, y = y, color = Income), shape = 17, size = 5) +
  coord_fixed()
ggplotly(p)

What is the mean income in zone 1? What is the mean income in zone 2? Not only are the summary measures of income highly representative of the observations they describe, the two zones are also highly distinct.

Imagine now that for whatever reason (lack of prior knowledge of the process, convenience for data collection, etc.) the zones instead are as follows:

zones2 <- data.frame(x1=c(0.2, 0.55), x2=c(0.55, 0.80), y1=c(0.0, 0.0), y2=c(0.6, 0.6), Zone_ID = c('1','2'))

If you plot these zones:

p <- ggplot() + 
  geom_rect(data = zones2, mapping = aes(xmin = x1, xmax = x2, ymin = y1, ymax = y2, fill = Zone_ID), alpha = 0.3) + 
  geom_point(data = df, aes(x = x, y = y, color = Income), shape = 17, size = 5) +
  coord_fixed()
ggplotly(p)

What is now the mean income of zone 1? What is the mean income of zone 2? The observations have not changed, and the generating spatial process remains the same. You will notice, however, that the summary measures for the two zones are more similar in this case than they were when the zones more closely captured the underlying process.

Visualizing area data: choropleth maps

The initial step when working with spatial area data, perhaps, is to visualize the data.

Commonly, area data are visualized by means of choropleth maps. A choropleth map is a map of the polygons that form the areas in the region, each colored in a way to represent the value of an underlying variable.

Lets use ggplot2 to create a choropleth map of population in Hamilton. Notice that the fill color for the polygons is given by cutting the values of POPULATION in five equal segments. In other words, the colors represent zones in the bottom 20% of population, zones in the next 20%, and so on, so that the darkest zones are those with populations so large as to be in the top 20% of the population distribution:

ggplot(Hamilton_CT) + geom_sf(aes(fill = cut_number(Hamilton_CT$POPULATION, 5)), color = NA, size = 0.1) +
  scale_fill_brewer(palette = "YlOrRd") +
  coord_sf() +
  labs(fill = "Population")

Inspecting the map above, would you say that the distribution of population is random, or not random? If not random, what do you think might be an underlying process for the distribution of population.

Often, creating a choropleth map using the absolute value of a variable can be somewhat misleading. As seen in the map above, the zones with the largest population are also usually large zones. Any process that you might think of will be confounded by the size of the zones. For this reason, it is often more informative when creating a choropleth map to use a variable that is a rate, for instance population divided by area to give population density:

pop_den.map <- ggplot(Hamilton_CT) + 
  geom_sf(aes(fill = cut_number(Hamilton_CT$POP_DENSITY, 5)), color = "white", size = 0.1) +
  scale_fill_brewer(palette = "YlOrRd") +
  labs(fill = "Pop Density")
pop_den.map

It can be seen now that the population density is higher in the more central parts of Hamilton, Burlington, Dundas, etc. Does the map look random? If not, what might be an underlying process that explains the variations in population density in a city like Hamilton?

Other times, it is appropriate to standardize instead of by area, by what might be called the population at risk. For instance, lets say that we wanted to explore the distribution of the population of older adults (say, 65 and older). In this case, normalizing not by area, but by the total population, would remove the “size” effect, giving a proportion:

ggplot(Hamilton_CT) + 
  geom_sf(aes(fill = cut_number((Hamilton_CT$AGE_65_TO_69 +
                                   Hamilton_CT$AGE_70_TO_74 +
                                   Hamilton_CT$AGE_75_TO_79 +
                                   Hamilton_CT$AGE_80_TO_84 +
                                   Hamilton_CT$AGE_MORE_85) / Hamilton_CT$POPULATION, 5)),
          color = NA, 
          size = 0.1) +
  scale_fill_brewer(palette = "YlOrRd") +
  labs(fill = "Prop Age 65+")

Do you notice a pattern in the distribution of seniors in the Hamilton, CMA?

There are a few things to keep in mind when creating choroplet maps.

First, what classification scheme to use, with how many classes, and what colors?

The examples above were all created using a classification scheme based on the quintiles of the distribution. As noted above, these are obtained by dividing the sample into 5 equal parts to give bottom 20%, etc., of observations. The quintiles are a particular form of a statistical measure known as quantiles, of which the median is value obtained when the sample is divided in two equal sized parts. Other classification schemes may include the mean, standard deviation, and so on.

In terms of how many classes to use, often there is little point in using more than six or seven classes, because the human eye cannot distinguish color differences at a much higher resolution.

The colors are a matter of style, but there are coloring schemes that are colorblind safe (see here).

Secondly, when the zoning system is irregular (as opposed to, say, a raster), large zones can easily become dominant. In effect, much detail in the maps above is lost for small zones, whereas large zones, especially if similarly colored, may mislead the eye as to their relative frequency.

Another mapping technique, the cartogram, is meant to reduce the issues with small-large zones.

Visualizing area data: cartogram

A cartogram is a map where the size of the zones is adjusted so that instead of being the land area, it is proportional to some other variable of interest.

Lets illustrate the idea behind the cartogram here.

In the maps above, the zones are faithful to their geographical properties. Unfortunately, this obscured the relevance of small zones. A cartogram can be weighted by another variable, say for instance, the population. In this way, the size of the zones will depend on the total population.

Cartograms are implemented in R in the package cartogram.

CT_pop_cartogram <- cartogram_cont(as(Hamilton_CT, "Spatial"), "POPULATION")

Mean size error for iteration 1: 5.93989832705797
Mean size error for iteration 2: 4.22003636161907
Mean size error for iteration 3: 3.22484467962666
Mean size error for iteration 4: 2.67349944741397
Mean size error for iteration 5: 2.39730255727434
Mean size error for iteration 6: 2.27353151789255
Mean size error for iteration 7: 2.21619358572193
Mean size error for iteration 8: 2.1621830846263
Mean size error for iteration 9: 2.10302324046727
Mean size error for iteration 10: 1.97693517836395
Mean size error for iteration 11: 1.83404151346793
Mean size error for iteration 12: 1.68895026525741
Mean size error for iteration 13: 1.48501500944163
Mean size error for iteration 14: 1.32015982874586
Mean size error for iteration 15: 1.25319335368508

Notice that the value of the function cartogram (i.e., its output) is a SpatialPolygonsDataFrame. This object needs to be converted to an object of class sf if we wish to use ggplot2 to visualize it:

CT_pop_cartogram.t <- st_as_sf(CT_pop_cartogram)

Plotting the cartogram:

ggplot(CT_pop_cartogram.t) + 
  geom_sf(aes(fill = cut_number(Hamilton_CT$POPULATION, 5)), color = "white", size = 0.1) +
  scale_fill_brewer(palette = "YlOrRd") +
  labs(fill = "Population")

Notice how the size of the zones has been adjusted.

The cartogram can be combined with coloring schemes, as in choropleth maps:

CT_popden_cartogram <- cartogram(as(Hamilton_CT, "Spatial"), weight = "POP_DENSITY")

Please use cartogram_cont() instead of cartogram().

Mean size error for iteration 1: 29.038428707007
Mean size error for iteration 2: 26.8741045865025
Mean size error for iteration 3: 25.1188420421205
Mean size error for iteration 4: 23.62794892945
Mean size error for iteration 5: 22.3098917140505
Mean size error for iteration 6: 21.106667145772
Mean size error for iteration 7: 19.9815506326308
Mean size error for iteration 8: 18.9106321110581
Mean size error for iteration 9: 17.8793922806811
Mean size error for iteration 10: 16.8783214422446
Mean size error for iteration 11: 15.9023828765327
Mean size error for iteration 12: 14.9493939574685
Mean size error for iteration 13: 14.0197712378282
Mean size error for iteration 14: 13.115878344782
Mean size error for iteration 15: 12.2410530253921

Tidy and restore the data:

CT_popden_cartogram.t <- st_as_sf(CT_popden_cartogram)

pop_den.cartogram <- ggplot(CT_popden_cartogram.t) + 
  geom_sf(aes(fill = cut_number(Hamilton_CT$POP_DENSITY, 5)),color = "white", size = 0.1) +
  scale_fill_brewer(palette = "YlOrRd") +
  labs(fill = "Pop Density")
pop_den.cartogram

By combining a cartogram with choropleth mapping, it becomes easier to appreciate the way high population density is concentrated in the central parts of Hamilton, Burlington, etc.

grid.arrange(pop_den.map, pop_den.cartogram, nrow = 2)

This concludes this chapter.